﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace MathLib
{
    public class Matrix
    {
        public double[,] m;
        private int rows;
        private int cols;

        public Matrix(int rows, int cols)
        {
            this.rows = rows;
            this.cols = cols;
            m = new double[rows, cols];

            set_to_zero();
        }

        public Matrix clone()
        {
            Matrix copy = new Matrix(Rows, Cols);

            for (int row = 0; row < Rows; ++row)
            {
                for (int col = 0; col < Cols; ++col)
                {
                    copy.m[row, col] = m[row,col];
                }
            }

            return copy;
        }

        public int Rows { get { return rows; } }
        public int Cols { get { return cols; } }

        public Vector get_row_vector(int index)
        {
            if (index < Rows)
            {
                Vector row = new Vector(Cols);

                for (int i = 0; i < Cols; ++i)
                {
                    row.v[i] = m[index, i];
                }

                return row;
            }
            else
            {
                throw new InvalidOperationException();
            }
        }

        public Vector get_col_vector(int index)
        {
            if (index < Cols)
            {
                Vector col = new Vector(Rows);

                for (int i = 0; i < Rows; ++i)
                {
                    col.v[i] = m[i, index];
                }

                return col;
            }
            else
            {
                throw new InvalidOperationException();
            }
        }

        public static Matrix operator +(Matrix a, Matrix b)
        {
            if (a.Rows == b.Rows && a.Cols == b.Cols)
            {
                Matrix c = new Matrix(a.Rows, a.Cols);

                for (int i = 0; i < c.Rows; ++i)
                {
                    for (int j = 0; j < c.Cols; ++j)
                    {
                        c.m[i, j] = a.m[i, j] + b.m[i, j];
                    }
                }

                return c;
            }
            else
            {
                throw new InvalidOperationException();
            }
        }

        public static Matrix operator -(Matrix a, Matrix b)
        {
            if (a.Rows == b.Rows && a.Cols == b.Cols)
            {
                Matrix c = new Matrix(a.Rows, a.Cols);

                for (int i = 0; i < c.Rows; ++i)
                {
                    for (int j = 0; j < c.Cols; ++j)
                    {
                        c.m[i, j] = a.m[i, j] - b.m[i, j];
                    }
                }

                return c;
            }
            else
            {
                throw new InvalidOperationException();
            }
        }

        public static Matrix operator*( Matrix A, Matrix B )
        {
            if (A.Cols == B.Rows)
            {
                Matrix c = new Matrix(A.Rows, B.Cols);

                for (int row_A = 0; row_A < A.Rows; ++row_A)
                {
                    for (int col_B = 0; col_B < B.Cols; ++col_B)
                    {
                        double inner_product = 0;

                        for (int element_index = 0; element_index < A.Cols; ++element_index)
                        {
                            inner_product += A.m[row_A, element_index] * B.m[element_index, col_B];
                        }

                        c.m[row_A, col_B] = inner_product;
                    }
                }

                return c;
            }
            else
            {
                throw new InvalidOperationException();
            }
        }

        public static Matrix operator *(Vector v, Matrix a)
        {
            if (v.N == a.Rows)
            {
                return v.get_matrix_as_col().tp() * a;
            }
            else
            {
                throw new InvalidOperationException();
            }
        }

        public static Matrix operator *(double s, Matrix a)
        {
            Matrix c = new Matrix(a.Rows,a.Cols);

            for (int row = 0; row < a.Rows; ++row)
            {
                for (int col = 0; col < a.Cols; ++col)
                {
                    c.m[row, col] = s * a.m[row, col];
                }
            }

            return c;
        }

        public void set_from_matrix(Matrix a)
        {
            if (Rows == a.Rows && Cols == a.Cols)
            {
                for (int row = 0; row < Rows; ++row)
                {
                    for (int col = 0; col < Cols; ++col)
                    {
                        m[row, col] = a.m[row, col];
                    }
                } 
            }
            else
            {
                throw new InvalidOperationException();
            }
        }

        public void set_to_zero()
        {
            for (int row = 0; row < Rows; ++row)
            {
                for (int col = 0; col < Cols; ++col)
                {
                    m[row,col] = 0.0;
                }
            } 
        }

        public void set_to_identity()
        {
            if (Rows == Cols)
            {
                for (int row = 0; row < Rows; ++row)
                {
                    for (int col = 0; col < Cols; ++col)
                    {
                        if (row == col)
                        {
                            m[row,col] = 1.0;
                        }
                        else
                        {
                            m[row,col] = 0.0;
                        }
                    }
                }
            }
            else
            {
                throw new InvalidOperationException();
            }
        }

        public Matrix tp()
        {
            Matrix transpose = new Matrix(Cols, Rows);

            for (int i = 0; i < Rows; ++i)
            {
                for (int j = 0; j < Cols; ++j)
                {
                    transpose.m[j,i] = m[i,j];
                }
            }

            return transpose;
        }

        private void GetMinor(Matrix src, Matrix dest, int row, int col, int order)
        {
            // indicate which col and row is being copied to dest
            int colCount = 0, rowCount = 0;

            for (int i = 0; i < order; i++)
            {
                if (i != row)
                {
                    colCount = 0;
                    for (int j = 0; j < order; j++)
                    {
                        // when j is not the element
                        if (j != col)
                        {
                            dest.m[rowCount,colCount] = src.m[i,j];
                            colCount++;
                        }
                    }
                    rowCount++;
                }
            }
        }

        private double CalcDeterminant(Matrix mat, int order)
        {
            // order must be >= 0
	        // stop the recursion when matrix is a single element
            if( order == 1 )
                return mat.m[0,0];
    
            // the determinant value
            double det = 0;
    
            // allocate the cofactor matrix
            Matrix minor = new Matrix( order - 1, order - 1);
    
            for(int i = 0; i < order; i++ )
            {
                // get minor of element (0,i)
                GetMinor( mat, minor, 0, i , order);
                // the recusion is here!
                det += Math.Pow( -1.0, i ) * mat.m[0,i] * CalcDeterminant( minor,order-1 );
            }
    
            return det;
        }

        public Matrix inverse()
        {
            if (Rows == Cols)
            {
                Matrix inv = new Matrix(Rows, Rows);
                inv.set_to_zero();

                if (Rows > 1)
                {
                    Matrix A = new Matrix(Rows, Rows);
                    A.set_from_matrix(this);

                    // get the determinant of A
                    double det = 1.0 / CalcDeterminant(A, Rows);

                    Matrix minor = new Matrix(Rows - 1, Rows - 1);
                    minor.set_to_zero();

                    for (int j = 0; j < Rows; ++j)
                    {
                        for (int i = 0; i < Rows; ++i)
                        {
                            // get the co-factor (matrix) of A(j,i)
                            GetMinor(A, minor, j, i, Rows);

                            inv.m[i, j] = det * CalcDeterminant(minor, Rows - 1);

                            if ((i + j) % 2 == 1)
                            {
                                inv.m[i, j] = -inv.m[i, j];
                            }
                        }
                    }
                }
                else
                {
                    inv.m[0, 0] = 1.0 / m[0, 0];
                }

                return inv;
            }
            else
            {
                throw new InvalidOperationException();
            }
        }
        
    }
}
